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^ ' Abstract 
I— 5 ! . 

, A physical model and two-dimensional numerical method for computing the evolution and 

I spectra of protostellar clouds are described. The physical model is based on a system of 

magneto-gasdynamical equations, including ohmic and ambipolar diffusion, and a scheme 
for calculating the thermal and ionization structure of a cloud. The dust and gas tem- 
peratures are determined during the calculations of the thermal structure of the cloud. 
^ ' The results of computing the dynamical and thermal structure of the cloud are used to 

Q^i model the radiative transfer in continuum and in molecular lines. We presented the results 



o 



for clouds in hydrostatic and thermal equilibrium. The evolution of a rotating magnetic 



■ protostellar cloud starting from a quasi-static state is also considered. Spectral maps for 
^ . optically thick lines of linear molecules are analyzed. We have shown that the influence of 

the magnetic field and rotation can lead to a redistribution of angular momentum in the 
cloud and the formation of a characteristic rotational velocity structure. As a result, the 
^ ■ distribution of the velocity centroid of the molecular lines can acquire an hourglass shape. 

'nI" , We plan to use the developed program package together with a model for the chemical 

^ I evolution to interpret and model observed starless and protostellar cores. 

O 1 Introduction 

a^ ■ 

^ ' Cool, dense, quasi-spherical condensations of gas and dust that do not contain infrared sources 
of radiation can often be distinguished in the observed complex structures of molecular clouds 
^ ' — so-called starless cores (see, for example [T]). These objects are characterized by densities of 

■ 10^-10^ cm~^, temperatures of 8-20 K, sizes of 0.1 pc, and masses from fractions of a solar mass 
to several tens of solar masses. It is believed that, in the absence of disruptive influences, they 
can enter a state of gravitational collapse that leads to the birth of young stars. In this case, 
these objects are usually called pre-stellar cores. Pre-stellar cores can evolve into protostellar 
cores, in which infrared sources are observed. We will use the term protostellar cloud as a 
gravitationally bound gas-dust condensation whose evolution leads to the formation of a single 
star or multiple star system. In various stages of its evolution, such a protostellar cloud can be 
manifest observationally as a starless or protostellar core. 

Starless and protostellar cores are observed primarily via the continuum (thermal) emission 
of dust and lines of compound molecules (CO, HCO"'", N2H"'" and others). Maps of dust emission 
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are used to establish the density distribution, while the molecular line emission provides infor- 
mation about the thermal, chemical, and dynamical structure of the cores (see, for example, 
[21 [I]). Both analytical (for example, the self-similar solutions [SHHE]) and numerical ([SIEIE]; 
see also the review [8]) models have been actively used in theoretical studies of the evolution of 
protostellar clouds. Numerical modeling has been used to study various effects and processes 
that influence the evolution of protostellar clouds, such as magnetic fields, turbulence, rotation, 
ambipolar diffusion, ionization, and chemical reactions. Moreover, numerical modeling is nec- 
essary for self-consistent computations of the characteristics of clouds based on observational 
data [9]. 

One promising direction is modeling the dust and molecular-line spectra and images of 
starless and protostellar cores based on modern dynamical models (see, for example, [TOl flT]). 
Since the line parameters (position, width, asymmetry) are sensitive to the velocity distribution 
inside the core, analysis of these lines can be used to determine the dynamical status of the core. 
At the same time, the molecular line emission depends strongly on many other parameters: the 
density distribution, temperature gradients, and concentrations of compound molecules [12]. 
This imposes additional requirements on the construction of dynamical models. A dynamical 
model for a starless or protostellar core should not only describe the evolution of the density 
and velocity distributions, but also adequately predict its thermal and chemical structure. 
So far, as a rule, modeling of the radiation spectra of protostellar clouds has been carried 
out using one-dimensional approximations. For instance, the evolution of the line profiles 
emitted by a protostellar cloud based on a one-dimensional dynamical model was studied in 
|13] . Theoretical molecular-line spectra for a chemical-dynamical model of the pre-stellar core 
L1544 are presented and compared with observations in [13]. 

Our aim here is to construct a physical model for starless cores and use it to model their 
spectra. Our approach includes a description of the main physical processes determining the 
dynamics of the core and the formation of the line profiles. Our dynamical model is two- 
dimensional and axially symmetrical, and takes into account the influence of rotation and 
magnetic fields, ohmic and ambipolar diffusion, ionization, heating and cooling processes, and 
radiative transfer. A non-LTE method is used to model the line radiation, which is able to 
describe both optically thin and optically thick line radiation. We emphasize that we rely 
on a complex approach that unifies a self-consistent dynamical model and a model for the 
radiative transfer in molecular lines. Such an approach has been used only in a one-dimensional 
approximation [TUl [HI [IS], which does not take into account important physical processes 
associated with rotation and magnetic fields. 

The article is organized as follows. Section 2 describes the main equations for the dynamical 
model and the numerical method used to solve them. Section 3 discusses the model for the 
radiative transfer in molecular lines that is used. Section 4 presents the results of numerical 
simulations of the structures of equilibrium and collapsing protostellar clouds. The Conclusion 
briefly discusses the main results of our study. 



2 



2 Dynamical model 



2.1 Basic equations 

We use the following system of equations to describe the dynamical evolution of magnetic, 
rotating protostellar clouds: 

|^ + V-(pv) = 0, (1) 
(9v 11 

— + (vV)v = — VP- — (B X (V X B)) - V<l>, (2) 
ot p Airp 

5B 

— = V X (v X B + vd X B - 7/oD (V X B)) , V • B = 0, (3) 



P ( § + (v ■ V) e ) + PV ■ V = - Ag - Agd, (4) 



= AnGp. (5) 

Here, p is the density, v the velocity, P the pressure, e the specific internal energy, B the 
magnetic field, and $ the gravitational potential. The induction equation (jS]) includes ohmic 
and ambipolar diffusion of the magnetic field. The ohmic diffusion coefficient is ?7od = c^/ (4vrcr), 
where a is the conductivity of the plasma. We have assumed standard magnetic ambipolar 
diffusion ([I5], Section 13.3), when the charged components of the plasma all move relative to 
the neutral components under the action of electromagnetic forces with essentially the same 
speed, 

(V X B) X B 

Vd = -r-B , 6 

where Ppn is the total coefficient of friction for the relative movement between the charged and 
neutral components. The energy equation (jlj) includes heating and cooling processes: Fg and 
Ag are the gas heating and cooling rates per unit volume and Agj the rate at which energy 
is redistributed between the gas and dust via collisions per unit volume. The system ([HIS]) is 
closed with the equation of state of an ideal gas, P = (7 — l)pe, where 7 = 5/3 is the adiabatic 
index. This value for the adiabatic index is dictated by the fact that the rotational degrees of 
freedom of the H2 molecules are frozen at the temperatures considered {5 K < T < 100 K) 



It is convenient to split the entire system ([UIH]) into simpler subsystems of equations accord- 
ing to the physical processes involved. The ideal magneto-gasdynamical equations can serve as 
the basic subsystem. Further, the Poisson equation for the gravitational potential, heating and 
cooling processes, radiative transfer, and ionization and diffusion of the magnetic field can be 
grouped into subsystems. Each subsystem can be solved using the most appropriate numerical 
methods. 

The dynamics of magnetic, rotating protostellar clouds can be studied in an axially sym- 
metrical approximation if the initial magnetic field B is colinear with the angular velocity of 
the rotation Q and there are no azimuthal perturbations. In this case, it is convenient to use 
cylindrical coordinates {r,ip,z). In the case of axial and equatorial symmetry, the numeri- 
cal modeling can be carried out in the two-dimensional computational domain {0 < r, z < R), 
where R is the size of the computational domain. The radius of the cloud was specified to be less 
than R, with the boundary of the cloud moving together with the gas during the compression 
process. 
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The conditions at the outer boundaries of the computational region r = R and z = R 
correspond to setting the derivatives of all computational quantities / in the direction of the 
outward normal to the boundary n equal to zero: df/dn = 0. The condition of equatorial 
symmetry is imposed in the z = plane, and the condition of axial sjTiimetry at the r = 
axis. 



2.2 Ionization balance equations 

To correctly describe the diffusion of the magnetic field in protostellar clouds, we must know 
the number densities of electrons rie, ions rii, and dust grains with charges of — e and +e (ng_ 
and rzg^, in accordance with [I?]). We used an ionization model based on that of [18], which 
is relatively simple, but provides good agreement with the results of complex astrochemical 
computations such as those of We included only one atomic ion, Mg"*", and one molecular 
ion, Hjj", with number densities of na+, Uj^^. The Mg"*" ion is formed in charge-transfer reactions 
between atomic Mg and Hj]", and also via photoionization by stellar radiation. The B.^ ion arises 
due to rapid reactions, beginning with the ionization of hydrogen and helium by cosmic rays, 
and the rate of production of B.^ is equal to the rate of ionization of hydrogen and helium, 
C = 3-10~^^ s~^ [20] • In the model used, the rates of ionization and recombination are sufficiently 
rapid that we can assume ionization equilibrium. Therefore, the equations for calculating na+, 
nra+, and n^^ have the form 

(-Pao + C'aom+'^m+) '^ao = (C'a+c^T-c + C'a+g.^. + C*a+go%o) ^g,^ , (7) 
C'^n = (C*dr^c + C*m+g_''^g_ + Cm+go^gg + C'aom+'^ao) ^ra+ y (8) 



Cego^cngo = [ Cg±ng^ + 2^C^+g_n^+ I ng_ , (9) 

'^go <^/3+go^/3+ = {Cg±ng_ + Ccg+rie) rig^ , (10) 
P 

where Pao is the probability of ionization of an Mg atom by stellar radiation per second (|15]. 
Section 5.2), riao = — n^^ the number density of neutral Mg atoms in the gas phase, ria 
the total number density of Mg in the gas phase, the number density of H, H2 and He, 
and rig^ the number density of neutral dust grains. The subscript (3 denotes Mg and H3. The 
notation for the reaction constants is indicated in the Table. The electron density is found 
from the condition of electrical neutrality. When calculating Pao? the photoionization cross 
section was set equal to the threshold value and taken out of the integral over frequency. The 
photorecombination coefficient Ca+e was taken to be that for a hydrogen-like ion ([ISj, Section 
5.1). We used the results of [22] to approximate the coefficient for dissociational recombination 
of H3 : Cdr = 1.28 • 10~^Tg~^''^ cm^s~^, where Tg is the gas temperature. 

In contrast to [18] we included the freezing of Mg onto dust grains. We assumed that, after 
the recombination of Mg"*" on a dust grain, the kinetic energy of the formed Mg atom exceeds 
the adsorption energy, so that only collisions with neutral atoms and dust grains were taken 
into account when calculating the Mg adsorption rate [2S]- At the gas and dust temperatures 
characteristic of protostellar clouds, the Mg desorption rate is much less than the Mg adsorption 
rate, so that the Mg abundance in the gas phase, Xa = ni^/n^, is appreciably non-equilibrium. 
Neglecting the convective derivatives in the continuity equations for the neutrals, Mg"*", and 
Mg, and also equating the divergences of their velocities, we obtain an approximate equation 
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Notation 


Reaction 


Reference 


C/3om+ 


Charge transfer for Hj]" and atoms (3 




C/3+e 


Photorecombination of ions (5 


m 




Recombination of ions (3 on negative grains 


m 


C/3+go 


Recombination of ions (3 on neutral grains 


m 




Dissociational recombination of 


[22] 




Electron capture by neutral grains 


HE] 




Mutual neutralization of grains 


dH] 




Neutralization of positive grains 


HE] 



Table 1: Notation and references for the ionization and recombination reaction constants 



that will be obeyed by Xa at a fixed point: 

d fn \ 

^ln^a= (^^-iJP.^a^gng (11) 

where Pg = 0.3 is the probability of adsorption of an atom during a collision with a dust grain 
(see, for example, [19]), fa the atom's mean thermal velocity, Sg and rig the mean geometrical 
cross section and number density of the dust grains. We used a high initial Mg abundance in 
the gas phase: Xa(0) = 5 • 10"'^. 

It is believed that the effective radii of dust grains have a power-law distribution with index 
—3.5 [15], with minimum and maximum radii of 10 nm and 1 /im. The ratio of the total 
concentrations of charged and neutral dust grains to Un is constant and determined from the 
mass fraction of dust, 1%. 



2.3 Description of the numerical method 
2.3.1 Magneto-gas dynamics 

The hyperbolic subsystem of equations separated out from the full system (IT]|5]) can be written 
in conservative form in cylindrical coordinates as follows: 



du_ 

'dt 



ld_ 

r dr 



dg 

dz 



(12) 



where the vectors of conservative variables U, fluxes JF and Q, and sources S are determined 
by the expressions 





p 






/ 


pVr 


\ 




pVr 








pvl + - p2/4vr 






pVz 








pVrVz — BrBz/4:7r 






pv^ 








pVrV^ - BrB^/A-K 






Bj. 













Bz 








VrBz — VzBr 






B^ 








VrB^ - V^Br 
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) 






:e + P,)-P,(v-B)/47r 


/ 



(13) 
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pv^Vr - B^Br/4:7r 
pvl + - BI/Ati 
pv^v^ - B^B^/An 

V^Br - Vr-Bz 



v^B^ - v^Bz 
V ^,(e + P,)-5,(vB)/47r / 



S 



( \ 

[pvl + - P^/47r) /r - pd^ldr 
-pd^/dz 
{-pVrV^ + BrB^/Air) /r 


(v^Br - VrBz) jr 


V -p(vV)$ ) 



(14) 



where P* = P + B^/Svr is the total pressure and, e = p (e + v^/2) + B^/Svr the total energy 
density. 

During the compression, there is an appreciable concentration of matter in the central part 
of the cloud, leading to strong radial gradients in the density, velocity, and other quantities. 
Therefore, to enhance the accuracy of the computations, we used an adaptive moving mesh 
in our numerical code, which became denser toward the center with time [21] . The system 
(I12til4p . was approximated in a moving curvilinear coordinate system that was related to the 
initial coordinate system by the transformations t = r, r = r(r, ^), z = z{t,()- The initial 
equations (11211141) in the new variables r, ^, and ( have the form 



o Q 

— {rQrQ^U) + — (rg,J^) + — [rQrG 



^QrQzS, 



(15) 



where T = T — uJrU-, G = G — ^zl^-, ^r,z is the velocity of the new coordinate system relative to 
the initial system, and Q^.z are metric coefficients. A more detailed description of the technique 
used to construct the adaptive mesh is presented in the Appendix. In the computational results 
presented below, the size of a central cell was approximately 0.1 AU, while the radius of the 
cloud itself was 0.2 pc. 

We introduced in the moving curvilinear coordinate system S,-, C in the computational domain 
a uniform difference mesh whose structure was determined by the distribution of the mesh 
coordinates = iA^, Q = jA^, where = R/Nj., A( = R/N^, and the subscripts i = 
0, 1, . . . , iVj. and j = 0, 1, . . . , A^^. The computational mesh quantities at time r" correspond 
to the cells c = (i + 1/2, j + 1/2), which are numbered with half-integer indexes. The values of 
the mesh quantities 

€2 C2 

1 



where Vr 



{rQrQ.YMAC 
(ri + r2)/2, at time r""^-*^ 



/ dC{rQrQMY 



(16) 



r 



«i Ci 

Ar were computed using the difference scheme 



At 



+ 



Ae 



AC 



{rQrQzS)l 



(17) 



The components of the source vector 5" were determined using expressions ( IHl) . which were 
computed at the center of cell c at time r". The operators Ag and A^ determine the difference 
between quantities at the cell boundaries in the corresponding directions. 

We computed the numerical vector fluxes JF and G through the cell boundaries by splitting 
them along spatial directions. The numerical fluxes in each spatial direction were computed 
based on the solution of the correspoding one-dimensional Riemann problem for the decay of 
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an arbitrary discontinuity. The split one-dimensional Riemann problem can be written in the 
form (we consider only the ^ direction as an example): 



Let us consider the solution of this problem in the approximation when all propagating waves 
that arise after the decay of the initial arbitrary discontinuity represent strong discontinuities. 
The corresponding Hugoniot conditions should be satisfied at each such discontinuity: 

DMa = [J" - cuMU (19) 

where a denotes the index of a discontinuity with velocity Da and the square brackets denote 
the difference between the right-hand and left-hand quantities for the given discontinuity. Note 
that, in this approximation for the solution of the Riemann problem, some propagating waves 
can represent non-evolving rarefaction shocks. 

The numerical fluxes JF through the cell boundaries in the scheme 0171) are determined as 
follows: 

^^ = - ^ E [^]" + (20) 

a 

where 

^^^ = IY1 1™^*^^" (^•^^-''i' ^^^^) - ^ E 1^^^^^^" (^•^^"' ^-^m) ' (21) 

a+ a— 

A^r^ = (22) 

The sum in the first term in (12T1) is carried out over all positive velocities Da, and the sum in 
the second term over all negative velocities. The function limiter(x, y) restricts the values of 
the anti-diffusion corrections SJ-' to ensure that the scheme is monotonic. We used the following 
non-linear limiting function for all a 



1 -\- p 1 — p 
limiter(x, y) = — - — minmod(i/, gx) H — minmod(gy, x), (23) 

where ^ 

minmod(x, y) = - [sign(x) + sign(|/)] min(|x|, \y\). (24) 

When p = 1/3 in regions where the solution is smooth, this scheme has third-order accuracy in 
the spatial variable and first-order accuracy in time. The value of q was equal to 4 in all our 
computations pU] . 

In our numerical computations, we used a difference scheme based on solving the Riemann 
problem in the HLLD approximation [27], which assumes that five strong discontinuities with 
four intermediate states are formed during the decay of the initial discontinuity: two MHD 
shocks, two Alfven (rotational) discontinuities, and an MHD contact (tangential) discontinuity. 
A generalized Lagrange multiplier (GLM) method was used to clean the divergence of the 
magnetic field [2H]- The described difference scheme for approximating the equations in the 
hyperbolic system (^S^ is a total variation diminishing (TVD) scheme with third-order accuracy 
in the spatial variable in smooth regions of the solution and first-order accuracy in time. With 
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the corresponding boundary conditions, the difference scheme approximating the ideal magneto- 
gasdynamical equations exactly conserves the total mass and energy. The stability of the 
scheme is provided by limiting the time step Ar (the Courant-Friedrichs-Levy condition) 
The Poisson equation for the gravitational potential was solved using ADI method [29] . 

To verify the computational properties and accuracy of the hyperbolic part of the numerical 
code, we carried out a series of test computations for problems with known exact or approx- 
imate analytical solutions (for details, see [221 EH])- Moreover, we checked the correctness of 
the operation of the code using self-similar solutions describing the inhomogeneous, free, and 
isothermal collapse of protostellar clouds [5]. The second self-similar solution corresponds to 
the critical case of the propagation of a rarefaction shock near the time when it is focused in the 
central part of the collapsing protostellar cloud. These computational results demonstrated the 
good computational properties of the code and its suitability for a wide range of astrophysical 
problems. 

2.3.2 Thermal Structure 

We separated out a subsystem of equations to describe the thermal structure of a protostellar 
cloud: 

de 

= Tg - Ag - Agd, (25) 
Td - Ad + Agd = 0, (26) 

where e is the specific internal energy of the gas, Fg and Fd the gas and dust heating rates per 
unit volume, and Ag and Ad the gas and dust cooling rates per unit volume. Equations fl25|) and 
( I26p are interrelated by the rate of the redistribution of energy between the gas and dust due 
to collisions between gas molecules and dust grains, Agd. Thus, the temperature distributions 
for the dust and gas must be found jointly with the solution of this system. The estimates of 
[30] show that, given the relatively low heat capacity of the dust, the time required to establish 
radiative equilibrium for the dust component is much shorter than the characteristic dynamical 
time-scale for a wide range of densities and temperatures in the models considered. Therefore, 
we neglected the time derivative of the thermal energy of the dust in 0261) : i.e., we assumed 
that the dust was in radiative equilibrium at each time. Further, we describe the heating and 
cooling processes that must be included when modeling the cores of molecular clouds. 

One of the main mechanisms for heating the gas, especially in the central, dense regions of 
the cloud, is heating by cosmic rays. We used the following formula from [31j : 

F„ = 10"2'n(H2), (27) 

where n{B.2) is the number density of molecular hydrogen. However, the intensity of the cosmic- 
ray flux is rather uncertain, and can vary appreciably for different clouds (see, for example, [32]). 
An important gas heating mechanism in the envelope is the photoelectric effect [331 when 
high-energy interstellar UV photons eject electrons from the surface of dust grains. These 
electrons have comparatively high energies, and so heat the gas. To calculate the heating 
function for the photoelectric effect, we used formula (42) from [S^ : 

Fpe = 10-2^peGon(H), (28) 

where epe is the efficiency of photoelectric heating, ?7,(H) the density of hydrogen (atomic and 
molecular), and Go the number of high-energy photons (6 eV < hu < 13.6 eV) at a given 
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point normalized to the mean number of high-energy photons in the interstellar radiation 
field (see [M], Section 4.1). 

We calculated epe using formula (43) from [33], taking the number density of free electrons 
in this formula from the ionization model described in Section 12.21 In our computations, epe 
varied from 0.05 at the center of the cloud to 0.01 in the cloud envelope. 

The gas in the cloud is primarily cooled by molecular-line radiation and collisions between 
the molecular gas and dust grains, which, as a rule, have lower temperatures. We used the 
following approximate formula from [33] to compute the molecular-line cooling function: 



where Tg is the gas temperature and the coefficients a and P depend on the density (see [35] . 
Table 2). This formula generalizes results of modeling the radiative transfer in a molecular 
cloud and of computing the cooling functions for cooling in lines of CO, ^^CO, C, and other 
compound molecules. The rate of energy exchange between the gas and dust was calculated 
using formula (15) from |35j : 

Agd = 6.32 ■ W-'^n\R2) T^^^ (Tg - Td), (30) 

where is the dust temperature. Depending on the ratio of and Tg, the energy exchange 
between the gas and dust can lead to either cooling and heating of the gas. 
The main heating mechanism for the dust is the absorption of radiation: 

oo 

Td = 47r J pdK^Judu, (31) 



where pd is the mass of dust per unit volume, k^, the absorption coefficient per unit mass of 
dust, and 

j^ = — [ i^dn (32) 



A-ir _ 

in 



the angle-averaged spectral intensity of the radiation. The dust is cooled primarily by its own 
blackbody (thermal) radiation: 

oo 

Ad = 47r j pd K^B^{Td) du, (33) 



where B^{T^) is the Planck function for the dust temperature Td. Computation of the mean 
intensity (l32ll requires knowledge of the spectral intensity, which can be obtained from the 
solution of the radiative transfer equation: 

^ = -pAi^vh + pdfiuB^{Td), (34) 
as 

where the first term on the right-hand side describes absorption of the radiation and the sec- 
ond term describes the blackbody radiation of the dust along a specified direction. Thus, 
determining the temperature of the dust requires a consideration of the system ( l26l) . ( 13T]) -(13 
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Figure 1: Scheme for integrating the radiative-transfer equation. The mean intensity is cal- 
culated using the short-characteristic method, where the rays along which the integration is 
carried out lie in the meridional plane. The integration in each direction is carried out to cell's 
boundary, where the intensity is found via interpolation. 

describing the condition for radiative equilibrium with an internal source, Agd. We now describe 
the method used to solve this system for the conditions characteristic of starless cores. 

First and foremost, we assume that the cloud is heated primarily by external (interstellar) 
radiation, and that the cloud itself is optically thin to its own (thermal) radiation; i.e., we 
neglect the second term in ( IMl) : 

-T" = -pdK,uIu- (35) 
as 

This equation describes the weakening of the interstellar radiation due to its propagation 
through the cloud. This approximation is well satisfied for molecular cloud cores, and sub- 
stantially simplifies the solution procedure. Further, we approximate the integral (l32l) in the 
form 

1 ^ 

1=1 

i.e., we consider a discrete set of N directions for which the intensities are determined. An 
important aspect of the organization of the computations is the choice of these directions. Since 
we are considering an axially symmetrical problem, we chose directions lying in the meridional 
plane (Fig. [T]). 

Of course, this choice introduces some small error, however it appreciably simplifies the 
solution procedure, since we need not consider curvilinear trajectories during the integration. 
We used the short characteristic method to solve the transfer equation for /i(z/) [32] • In this 
method, the integration (l35l) in each cell is carried out from its center to its boundary, with the 
intensity at the cell boundary found via interpolation. To model the radiative transfer, we used 
the frequency dependence of the opacity n^, for compact, silicate dust grains with ice mantles 
from [37] . 

We tested the described method for modeling the thermal structure of the cloud using 
a series of model problems. In particular, we compared our computational results with the 
solution obtained by the TRANSPHERE-llf] numerical code developed by C. DuUemond to 

"'^http: / / www.mpia.de /homes/dullemon/radtrans 
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0.001 0.01 0.1 0.001 0.01 0.1 



r, pc r, pc 

Figure 2: The ratio of the heating/coohng rates of the gas to the heating rate by cosmic rays 
for equihbrium models of protostellar clouds with central densities Uc = 10^ cm"^ (left) and 
Uc = 10^ cm~^ (right) 

model radiative transfer in dust envelopes taking into account the thermal radiation of the dust. 
The maximum difference in the dust temperatures in the obtained solutions was less than 0.2 K 
for all conditions considered, justifying our use of the approximation 

To demonstrate the relative contributions of various processes to the thermal balance of the 
cloud, we present the ratio of heating and cooling rates of the gas to the heating function by 
cosmic rays for the equilibrium cloud model described in detail in Section 14.11 Figure [2] shows 
our results for models with central densities ric = 10^ cm~^ and = 10^ cm~'^. 

In the model with Uc = 10^ cm~^, the heating by cosmic rays and the photoelectric effect 
are comparable in magnitude in the central regions of the cloud, while photoelectric heating 
dominates in the cloud envelope. Cooling of the gas due to collisions with dust is appreciably 
lower than cooling via molecular-line radiation in the entire volume of the cloud. The dominant 
heating mechanisms in the cloud center and envelope in the model with Uc = 10^ cm~^ are 
clearly distinguished: cosmic rays and the photoelectric effect, respectively. The cooling mech- 
anisms are also clearly distinguished: collisions with dust in the cloud center and molecular 
line radiation in the envelope. 

2.3.3 Diffusion of the magnetic field. 

To describe the diffusion of the magnetic fields in collapsing protostellar clouds, we can neglect 
the friction between the charged components of the plasma [17]. We therefore assume that a 
component of type a, consisting of particles with mass ttIq,, charge ±e, and number density Ua 
has conductivity cTq, = e'^n^/Ran, where Ran = fJ'annann{sv) an is the coefficient of friction with 
neutral particles, ^an the reduced mass, and {sv)an the deceleration coefficient. Dividing the 
charged particles into four types — electrons, ions, and positive and negative dust particles 
— we can write the total conductivity of the plasma a = ac + di + cig^ -|- (Tg_ , and the total 
coefficient of friction between the charged and neutral particles R^n = -Ren + Rin + -Rg+n + -Rg_n- 
For the ions, {sv)in = 2 ■ 10"^ cm^ ■ s~^ ([IS], Section 2.2). When computing the deceleration 
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coefficient for the electrons, (sf)en we used the approximate tabulated data of ([38], p. 393). 
We assume for the charged dust grains (s'y)g^n = {sv)g_n = SgV^ [IE], where is the mean 
thermal velocity of the neutral components. 

Using ([6]) the diffusion part of the induction equation ([3]) can be rewritten: 

— = r]V^B + V X (uadB) - V?] x (V x B), (37) 

where r/ = ?7oD + B^/(47ri?pn) is the total diffusion coefficient and mad = (B ■ (V x B))/(47r_Rpn) 
the scalar ambipolar diffusion velocity. Since the diffusion coefficient rj and the scalar ambipolar 
diffusion velocity mad depend on the magnetic field, the use of explicit schemes to solve (1371) 
leads to strong restrictions on the time step. Therefore, to approximate 0371) in the interval 
r"" < r < r""*"^ in the two-dimensional numerical code, we used an absolutely stable, fully 
implicit scheme, which was realized via the following iterative process: 

r(p+i) _ r(o) . , , . 

= ?7(p)v'b(p+^) + V X (wilB^P)] - Vr/(P) x (V x B^^)), (38) 

where p is the iteration number. The values obtained via the solution of the hyperbolic sub- 
system of equations are used as the initial conditions for B*^^''. A splitting according to spatial 
directions (a locally one-dimensional scheme) is carried out in ( l38l) . and all differential operators 
are approximated with finite differences with second order accuracy. The resulting system of 
linear algebraic equations was solved with Thomas algorithm along each spatial direction. 

The ionization balance equations fITl fTIIl) were solved numerically using the Newton method. 
The resulting system of algebraic equations was solved using the method of Gauss. 



3 Model for radiative transfer in molecular lines 

As a rule, the formation of line radiation in molecular cloud cores occurs under non-LTE 
conditions (see, for example, [12] )• Therefore, modeling these lines requires the consideration 
of a system of balance equations for the level populations of the molecules, together with the 
radiative-transfer equation. We solved this problem using the approach described in detail in 
[9]. Here, we present only the main ideas of this method. 

The modeling of line radiation includes two stages. First, the distribution of the level 
populations (excitation temperatures of the transitions) of the molecules considered consistent 
with the external radiation and thermal radiation of the cloud was found. Accelerated A 
iterations were used for it, where the mean intensity in each cell of the computational domain 
was found using the long-characteristic method. The idea behind the accelerated A iterations is 
to use the radiation field computed in the previous iteration step to find a new approximation 
for the excitation temperature, which, in turn, is used to compute the intensity in the following 
iteration. The rays along which the integration of the transfer equation is carried out were 
chosen using the Monte-Carlo method. The integration of the radiative-transfer equation took 
into account the continuum radiation of the dust, with the dust temperature taken to be 
specified. 

In the second stage of the modeling, the obtained excitation-temperature distribution for 
the transitions was used to compute the observed properties of the modeled object, namely 
the line profiles, spectral maps, integrated intensity maps, mean velocity maps, and so forth. 
The orientation of the modeled cloud in space can be chosen, as well as the parameters of the 
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radio-telescope beam. Our numerical code also enables a detailed analysis of the formation of 
a studied line. The modeling results were used to calculate the contribution function which 
provides a relative input of all elements of the cloud to the resulting line profile [12]. Such 
analyses represent a powerful tool for interpreting the modeling results. 

The radiative transfer in molecular lines was modeled using a dedicated spatial grid that 
differs from that used to compute the dynamical evolution of the cloud. This is due to the fact 
that the criteria for the radiative transfer grid differ from those required for the hydrodynamical 
computations. Therefore, to construct the spatial grid for the line modeling we interpolated 
values from the hydrodynamical grid. The molecular data, such as the energy levels, Ein- 
stein coefficients, and collisional-excitation coefficients, were taken from the Leiden LAMBDA 
databas^ [39]. The frequency dependence of the opacity used for the dust radiation was the 
same as that used when modeling the thermal structure of the cloud. The code was well tested, 
and successfully used in the past to model radiative transfer in molecular cores [101 [12] and 
protoplanetary disks [HI Il2] . 

4 Demonstrative computations 

4.1 Equilibrium configurations of protostellar clouds 

The equilibrium structure of a protostellar cloud can be described using the system of equations: 



where TZ is the universal gas constant and the mean molecular weight. These equations were 
solved numerically in a spherically symmetrical approximation, iterating in the temperature of 
the gas and dust until the required accuracy was attained. The equations of hydrostatics (l39l) 
and thermal balance of the gas and dust fHOj) were solved in each iteration. The thermal-balance 
equation was solved using the Newton method, jointly with the radiative-transfer equation 
( I35p and the ionization-balance equations flTl fTUl) . The densities at the cloud center and in the 
outer envelope were specified as boundary conditions. The radius and mass of the cloud were 
determined from the solution of ( l39l l40l) . 

Figure [3] shows the distributions of the number density of molecular hydrogen n(H2), the 
dust temperature, and the gas temperature for equilibrium cloud models with central densities 
Tie = 10^ cm~^ (left) and = 10^ cm~^ (right). The envelope density is riext = 10'^ cm~^ in 
both models. 

The radius and mass of the cloud core in the first model are -Rcore = 0.2 pc and Mcore = 
4.4 Mq. The temperatures of the gas and dust at the cloud center are Tg = 12.3 K, Td = 11 K, 
while the temperatures in the outer envelope are Tg = 33.1 K, = 16.9 K. The difference in 
the gas and dust temperatures in the central parts of the cloud is due to the cooling of the gas 
in the cloud via molecular line radiation (left panel in Fig. |2]). 

In the second model, -Rcore = 0.2 pc and Mcore = 4.6 Mq. The gas and dust temperatures 
at the cloud center are Tg = 5 K and Td = 4.1 K, while the temperatures in the envelope are 
Tg = 33.5 K and Td = 17 K. The gas and dust temperatures in the central regions of the cloud 
are approximately the same, since the gas is mainly cooled via radiation by dust (right panel in 

■^http: / /www. strw.leidenimiv.nl/~nioldata/ 



VP = -pV$, V2<l> = 4nGp, P = 7^pTg//i, 



(39) 



Fg - Ag - Agd = 0, Fd - Ad + Agd = 0, 



(40) 
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Figure 3: Distributions of the number density of molecular hydrogen, the gas temperature, and 
the dust temperature for equilibrium cloud models with central densities Uc = IC^ cm~^ (left) 
and ric = 10^ cm~^ (right). 

Fig. [2]). In the second model, a central, cool (with a mean gas temperature of about 7 K) cloud 
core with a radius of about 0.05 pc and a mass of 1 Mq and an extended thermal envelope with 
a characteristic gas temperature of 12 K are clearly distinguished. 

The cloud masses and radii are the same in the two models. This means that these models 
can be taken to be different stages in the quasi-static evolution of a single cloud. The equilibrium 
configurations obtained can be used to analyze the structures of observed starless cores, since 
they are more accurate than the widely used isothermal Bonnor-Ebert models. 

4.2 Collapsing protostellar clouds 

As an example, we describe the results of our numerical modeling of the compression of a 
protostellar cloud with initial parameters corresponding to the first equilibrium configuration 
from the previous section (central density Uc = 10^ cm""^, envelope density riext = 10^ cm~^, 
radius -Rcorc = 0.2 pc, mass Mcove = 4.4 Mq). We specified in addition a uniform magnetic field 
{Bq = 2.6 fiG, ratio of magnetic energy to the modulus of the gravitational energy 0.05) and 
rigid rotation (Oq = 3.3-10"^^ s~^, ratio of rotational energy to the modulus of the gravitational 
energy 0.1) in the cloud. This uniform magnetic field is force-free, and so does not influence 
the equilibrium structure. The rotation creates a centrifugal force, which perturbs the initial 
equilibrium state. However, we specified fairly slow rotation of the cloud, which did not lead to 
appreciable deviations from the equilibrium state. As the computational results showed, this 
demonstration model has rich physical properties that are clearly manifest in the molecular line 
profiles. 

Figure H] shows the distributions of the logarithm of the density (gray scale and contours) 
and the poloidal velocity (arrows) in the cloud at time t = 5.8 tff = 2 million years from the 
initial contraction, where t// = l/y/AirGpo is the free-fall time and, po the initial density at the 
cloud center. By this time, the central density has reached ric = 10^ cm~^. The velocity of the 
collapse is slowed near the equatorial plane by the action of the electromagnetic and centrifugal 
forces. Therefore, the cloud takes on an elongated structure with time. 

Figure E] shows distributions of the abundances of charged particles in the equatorial plane 
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Figure 4: Modeling results for the collapse of a protostellar cloud. The left panel shows the 
distributions of the logarithm of the density (gray scale) and the poloidal velocity (arrows), 
while the right panel shows the distribution of the rotational velocity (gray scale) and the 
magnetic field lines. 



at t = and t = 5.8 tjf. In both cases, Mg"*" is the dominant ion in the outer layers of the 
cloud, where the photoionization of Mg atoms is substantial and the adsorption of Mg atoms 
onto dust grains is weak. In the dense central layers of the cloud, adsorption of Mg occurs on 
a time scale comparable to or shorter than tff, and the abundance of Mg"*" becomes negligibly 
small by t = 5.8 tff. The abundance of in our model is appreciably higher in all layers than 
in more complex models, such as that of [I9], providing a good agreement with the electron 
abundance. The abundances of charged dust grains remain nearly constant in the collapse 
stages considered, and are appreciably lower than the abundances of ions and electrons. 

Magnetic ambipolar diffusion leads to appreciable damping of MHD waves, but only weakly 
influences the evolution of the large-scale magnetic field. Therefore, a toroidal magnetic-field 
component, is generated in the collapsing cloud due to its differential rotation. This gives 
rise to a braking torque that brings about a redistribution of angular momentum between the 
central and outer regions of the cloud. The right panel in Fig. H] shows the distribution of 
the azimuthal velocity (gray scale) and the magnetic field lines at time 5.8 tff. This shows 
that the region of the redistribution of angular momentum is determined by the magnetic tube 
formed by the toroidal magnetic field. The local specific angular momentum is redistributed 
within the cloud, among its central and peripheral regions, and is also carried into the external 
medium along the formed magnetic tube. The opening angles of the magnetic tube and the cone 
of the maximum azimuthal velocity nearly coincide, and are approximately 40°, in agreement 
with theoretical estimates and observed opening angles for magnetic field lines in starless and 
protostellar cores [3]. 

Figure [6] presents the results of modeling the radiative transfer in the HCO'^(l-O) line, 
which is often considered in observations of starless cores. We assume that the number density 
of HCO''' relative to hydrogen is constant throughout the cloud and equal to 10~^, which 
corresponds to the typical mean values in starless cores [12]. We consider the case when the 
inclination of the symmetry axis of the cloud relative to the observer is i=60°, and the rotation 
of the projected axis in the plane of the sky is PA=0°. The distance to the cloud is taken to 
be 140 pc. We have not convolved the results with any telescope beam. We can see on the 
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Figure 5: Distributions of the abundances of charged particles in the equatorial plane for a 
model of a collapsing cloud for t = 0, ric = 10^ cm~^ (left) and t = 5.8 tff, Uc = 10^ cm~^ 
(right) 



spectral map (left) that the line profiles have a complex, asymmetric shape with this asymmetry 
changing with position on the plane of the sky. The formation of these asymmetric profiles is 
due to the self-absorption of radiation in an optically thick medium with a velocity gradient 
along the line of sight [12] . Profiles with double peaks dominate in the right-hand part of the 
map, with the left peak stronger than the right peak. The right peak is barely visible in the 
line profiles in the left-hand part of the map. This can be understood if the projections of the 
infall and rotation velocities along the line of sight have different signs in the right-hand part 
of the map, leading to the formation of more symmetrical, two peaked profiles, while these two 
projected velocities have the same sign in the left-hand part of the map, so that the right peaks 
are essentially fully suppressed. The map of the integrated intensity (right) is also asymmetric. 
In particular, the position of the intensity maximum does not coincide with the direction toward 
the center of the cloud, due to the described variation in the line asymmetry. 

Figure [7] presents maps of the velocity centroid (the first moment of the line profile or 
mean velocity) for the central regions of the cloud for i=QO° (left) and for the entire cloud for 
i=90° (right). A gradient in the mean velocity due to the rotation of the central regions of the 
protostellar cloud is clearly visible in the left panel. At the same time, the mean velocity is 
predominantly negative, testifying to a collapse, in accordance with the velocity distributions 
in the model. The rotational-velocity field has an hourglass shape (right panel in Fig. Hj) 
and is clearly manifest in the distribution of the velocity centroid for the entire cloud (right 
panel in Fig. [7]). At large distances from the cloud center, where the manifestations of the 
collapse are weaker, the distribution of the velocity centroid becomes more symmetrical about 
the projection of the rotational axis (white and black contours). Like the opening angle for 
the cone representing the maximum azimuthal velocity in the dynamical model (right panel in 
Fig. HI) the opening angle for the velocity contours far from the cloud center is roughly 40°. The 
described observational manifestations of the rotation of the cloud can be used when planning 
observations of starless cores. If these properties are detected, this will provide evidence for the 
mechanism for carrying away angular momentum described above. 
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Figure 6: Results of modeling the radiative transfer in the HCO'^(l-O) line for a collapsing 
protostellar cloud. The left panel shows a map of the line profiles for the central regions of the 
cloud, while the right panel shows a map of the integrated radiation intensity in K km/s. The 
distance to the cloud is 140 pc, the inclination of the symmetry axis relative to the observer 
i=60°, and the rotation of the projected axis in the plane of the sky PA=0°. The asymmetry 
of the line profiles is due to the collapse velocity and the differential rotation of the cloud. 



5 Conclusion 

We have presented a physical model and corresponding program package for the computation 
of the evolution of protostellar clouds and their observational properties. The physical model 
is based on magneto-gas-dynamical equations and methods for computing the thermal and 
ionization structure of a cloud. Ohmic and ambipolar diffusion are taken into account, as well 
as radiative transfer through dust. The calculated dynamical and thermal structures of clouds 
are used to model the radiative transfer in molecular lines. We have presented the results of 
modeling the structure of a quasi-equilibrium protostellar cloud and its subsequent evolution 
as an illustrative example. 

We have not considered the chemical structure of starless cores, which could substantially 
influence the formation of the molecular-line profiles. Analysis of observations of starless cores 
indicates that many molecules are frozen onto the surfaces of dust grains in the central regions 
[121111], while molecules are destroyed by interstellar UV radiation in the cloud envelopes [ID]. 
Therefore, it is important to use a self-consistent chemical-dynamical model for modeling the 
corresponding line profiles and interpreting observations. The development of a module for 
computing the chemical evolution is a separate task that will be the subject of future work. 
After the addition of this chemical module, it will be possible to use our program package not 
only for theoretical studies of the structure of starless cores, but also for modeling of particular 
observed objects. 

Our program package can also be used to model more advanced stages of the collapse of 
protostellar clouds, right up to the formation of an opaque core (protostar). However, this 
requires modification of the method used to compute the radiative transfer for the case of an 
optically thick medium. 



17 




150 100 50 -50 -100 -150 400 200 -200 -400 

RA, orcsec RA, orcsec 



Figure 7: Results of modeling the radiative transfer in the HCO+(1-0) line for a collapsing 
protostellar cloud. The left panel shows a map of the velocity centroid for the central regions 
of the cloud (the inclination angle is i=60°). while the right panel shows a map of the velocity 
centroid for the entire cloud (z=90°). The velocity is given in m/s. 
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Appendix 

A Adaptive mesh 

Let us describe the method used to construct the adaptive, moving mesh used in the two- 
dimensional numerical code. The theory of such meshes is based on canonical transformations 
of a hyperbolic system of equations in conservative form. As applied to the system fll2m4p . we 
will use the following transformation of coordinates and time: 

t = T, r = r(r,0, z = z{tX), (41) 
which can be rewritten in differential form 

dr = Uj-dr + QrdC,, dz = uj^dr + QzdC, (42) 
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where ujr,z is the velocity of the mesh and Q^.z are metric coefficients. The law for the trans- 
formation of the coordinates fl41l) is chosen in the form 

t = r, r(r, = ^ g^^^^n/A^ _ ^(^' C) = R^j^yi/K^- (^3) 

Using these relations, we can obtain explicit expressions for the velocities of the mesh Ur and 
(jUz and the metric coefficients Qr and Q^. 

The positions of the nodes rj(r), Zj^r) depend on time, so that, generally speaking, the 
mesh in the initial coordinates (r, z) at an arbitrary time is inhomogeneous. It is not difficult 
to convince oneself that the transformation law ( H3l) corresponds to the conditions 



where Ari+i/2 = rj+i — and Azj+1/2 = zjj^i — zj are the mesh steps in the r and z directions. 
Thus, the quantities qr,z{'T) > 1 are parameters of the mesh compression in the corresponding 
directions. When q^^z = 1? the mesh in (r, z) is uniform, since then rj = ^j, Zj = Q. 

The time dependence of the compression parameter can be specified arbitrarily taking into 
account the assumed character of the flow. The contraction of protostellar clouds has a char- 
acteristic dynamical time scale — the free-fall time t// = 1 / y/^TrUpo, where po is the initial 
density at the cloud center. We used the following expression for the compression parameter: 

q{T) = goo - (goo - go) (l + r/r,) (45) 

where go, goo, and r* are parameters. This dependence was chosen so that g(0) = qo, q ^ goo as 
r — > oo, and variations of the function g(r) occur on the characteristic time scale r*. Moreover, 
g(0) = 0. The initial and final values of the compression parameters qr^r) and qzi^) should 
not be very different. We chose these so that the relative difference in the sizes of neighboring 
mesh cells in (r, z) did not exceed 10-20%. The characteristic time scale was chosen to be 
several free-fall times tff. 
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